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The  RMS/MTTS  Instrumentation  system  located  at  MacGregor  Range  is  a  range 
measuring,  multiple  target  tracking  system.  In  order  to  obtain  a  vehicle 
trajectory  from  this  system,  the  range  measurement  from  several  receivers  are 
processed  by  least  squares.  Since  the  measured  vehicle  trajectories  are  often 
low  altitude,  the  resulting  nonlinear  least  squares  equations  are  111 -condi tione< 
In  addition,  this  measurement  system  Is  plagued  by  outliers,  sometimes  by  dense 
burst  of  outliers.  The  combination  of  ill -conditioning  and  outliers  is  lethal 
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and  attempts  to  robustify  the  nonlinear  least  squares  processing  have  failed. 
An  alternative  is  to  preprocess  each  of  the  range  measurement  sequences, 
eliminating  the  outliers  and  replacing  them  if  necessary.  Each  sequence  of 
range  measurements  is  preprocessed  by  robustly  fitting  a  cubic  spline 
using  iteratively  reweighted  least  squares.  Due  to  the  nature  of  spline 
fitting  and  the  possible  dense  bursts  of  outliers,  the  choice  of  a  good  set 
of  initial  weights  for  use  in  the  iteratively  reweighted  least  squares  is 
important  to  the  efficiency  of  the  method.  These  initial  weights  are 
determined  using  robust,  local  fitting  techniques.  Several  robust  techniques 
have  been  tested  for  this  local  fitting  application.  The  robust  spline 
preprocessing  Is  illustrated  with  some  especially  troublesome  data  sequences 
and  the  relative  performance  of  several  robust  methods  for  choosing  the 
initial  weights  is  compared.^ _ 
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The  RMS/MTTS  instrumentation  system  located  at  MacGregor  Range  is  »  range 
measuring,  multiple  target  tracking  system.  In  order  to  obtain  a  vehicle 
trajectory  from  this  system,  the  range  measurements  from  several  receivers 
are  processed  by  least  squares.  Because  the  measured  vehicle  trajectories 
of  Interest  are  often  low  altitude  and  because  of  the  geometry  of  the  re¬ 
ceiving  stations,  the  resulting  nonlinear  least  squares  equations  are  o*ten 
ill-conditioned.  In  addition,  this  measurement  system  Is  subject  to  out¬ 
liers,  sometimes  dense  bursts  o *  outliers.  This  combination  of  ill-condition¬ 
ing  and  outliers  Is  lethal  and  our  attempts  to  robustlfy  the  nonlinear  least 
squares  estimation  process  have  failed.  An  alternative  Is  to  preprocess  each 
of  the  range  measurement  sequences.  Identifying  the  outliers,  and  replacing 
them  if  necessary.  The  Ill-conditioned  least  squares  problem  can  then  be 
treated  without  being  troubled  by  outliers. 


Suppose  we  preprocess  the  measurement  sequence,  R(t^),  1  *  1,N.  For  typical 
aircraft  trajectories  the  measurement  rate  is  10/sec  with  a  time  of  interest 
of  40  -  120  sec  so  that  N  Is  often  in  the  range  400  -  1200.  The  purpose  of 
the  preprocessing  may  be  to  detect  outliers,  to  precompute  measurement  variances 
for  future  least  squares  processing,  or  to  synchronize  several  different  dis¬ 
crete  measurement  sequences.  The  preprocessing  of  the  range  measurement 
sequence,  R(t^},  Is  done  by  fitting  a  cubic  spline  to  the  discrete  measurements 
using  Iteratively  reweighted  least  squares  (IRWLS).  Specifically,  at  the  k— 
Iteration  we  minimize. 


(1) 


1 


where  B. (•)  are  the  cubic  3-sp1ine$  and  bi  '  are  the  spline  coefficients  to 
be  estimated.  The  weights,  are  computed  from  the  Hampel  ^-function 
using  the  spline  fit  from  the  (k-1)—  iteration. 


where 


|>(x) 


/  x  |x|sa 
a*sgn  (x)a<|x|ib 

.(x-ch;r(<1) 


(3) 


s5k_1)  estimates  the  dispersion  in  the  residual,  R(t.)  -  J  b(k_1 V (t .) .  The 

«  .  j  i  i  i  j 

(k-1)  1 

value  of  Sj  '  can  be  computed  either  locally  or  globally  from  the  residuals 

at  the  (k-1 )— iteration.  The  dispersion  is  a  MAO  estimate  obtained 

from 

s|k_1)  -  median  |R(tm)  -  J  b\M  \  (tm)  |/.6745  (4) 

J  meT.  i  1 

If  the  set  Tj  is  in  some  sense  the  set  of  points  close  to  tj,  the  estimate 
1S  1ocal*  If  set  Tj  the  set»  Tj  *  {tjm  »  1  ,N>  the  estimate 

is  global.  For  the  present  application  only  the  global  estimate  S^k_^«  s^k’^ 
will  be  used.  If  a  very  long  data  sequence,  say  about  one  hour,  a  local 
estimate  would  probably  be  preferable  to  the  global  estimate. 

CHOICE  OF  KNOTS 

Let  {T^ ,  i  »  1 ,M}  be  a  set  of  knot  times.  These  knot  times  are  used  to  define 
the  cubic  8-spllnes,  Bj(t^).  Of  most  importance  In  the  choice  of  the  knot  times 
is  their  spacing,  which  determines  the  ability  of  the  cubic  spline  to  fit  the 
data.  However,  for  each  additional  knot  time  there  Is  one  additional  spline  co¬ 
efficient  to  be  estimated,  thus  Increasing  the  computational  load.  Thus,  we  want 

• 

to  have  as  few  knots  as  possible  and  the  rules  for  their  choice  simple  and  yet  be 
able  to  adequately  represent  the  data.  With  this  simple  philosophy  for  selecting 
knots  we  will  try  to  assign  a  fixed  number  of  data  points,  NPTO,  to  each  knot  in¬ 
terval.  The  first  four  knots  are  placed  at  the  first  data.  The  last  knot  inter¬ 
val  may  have  more  than  NPTO  points  but  fewer  than  2  •  NPTO  data  points.  If  there 
is  a  large  time  break  In  the  data,  a  knot  is  placed  at  the  beginning  and  end 


of  the  time  break.  The  Interval  between  these  two  knots  has  zero  data  oolnts 
and  the  interval  immediately  preceding  the  time  break  may  have  more  than 
NPTO  points  but  fewer  than  2  •  NPTO  points.  If  immediately  after  a  time  break 
there  is  a  second  time  break  before  NPTO  points  have  been  read,  the  few  (less 
than  NPTO)  points  read  between  the  two  time  breaks  are  discarded.  If  a  time 
break  occurs  while  reading  points  for  the  first  knot  interval,  the  few  (less 
than  NPTO)  points  are  discarded  and  the  first  four  knots  repositioned  at  the 
first  data  time  after  the  time  break.  If  a  time  break  occurs  during  the  last 
interval,  the  portion  of  the  last  interval  contiguous  to  the  previous  interval 
Is  kept  and  the  remainder  of  the  points  in  the  last  interval  are  discarded. 

If  there  are  at  least  NPTO  points  kept,  these  points  form  the  last  interval. 

If  there  are  less  than  NPTO  points  kept,  these  points  are  appended  to  the  pre¬ 
vious  Interval  so  that  the  number  of  knot  intervals  is  reduced  by  one.  The 
time  difference  between  successive  data  points  which  Is  used  to  define  a  time 
break  is  named  FITBRK.  FITBRK  is  dependent  on  the  sample  rate.  The  time  dif¬ 
ference  between  successive  data  points  used  to  define  a  time  break  in  the  first 
and  last  knot  intervals  is  FITBRK/5.  This  smaller  value  is  used  in  the  first 
and  last  Interval  because  It  is  critical  to  obtain  a  good  fit  in  these  intervals. 
The  flow  chart  on  the  following  pages  more  clearly  defines  the  logic  for  select¬ 
ing  the  knot  times.  The  following  define  the  variables  in  the  flow  chart: 


NOTS 

■  number  of  interior  knots 

R(-) 

*  array  of  range  measurements 

KR 

*  number  of  knot  Intervals 

NPTS(-) 

*  array  of  point  counts  for 

knot  intervals 

TT(  •) 

■  array  of  knot  times 

IBRK 

■  logical  denoting  the  occurrence 

NC0UN 

■  total  data  point  count 

of  a  time  break 

T(  •) 

■  array  of  data  times 

STA 

*  data  start  time 

DO  FOR  I  »  1 ,  N 


KR  -  NOTS  -  1 
! INDX  «  NCOUN  -  NPTS(KR) 
TL  -  T( I NDX  +  I ) 


DO  FOR  K  -  2,  NPTS(KR) 


INDX  +  K)  >  .2*FITBR 


<  NPTO>— Sj NOTS  »  NOTS-1  kwS  ^V>  THEN 


'ELSE 


NPTS(KR)  »  K-l 


KR  *  KR  -  1 
N?TS(KR)  »  Vots(<(R)+K 


TT(NOTS  +3)  *  TL  +  EPS 
ETA  -  TL  _ 


TL  *  T(INDX)  +  K 
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THE  LEAST  SQUARES  NORMAL  EQUATIONS 

At  the  k—  iteration  of  the  fitting  procedure  the  weighted  sum  of  squares 


I  uik)  (R(t,)  -  I  b,(k)  B,(t,))‘ 
j*l  J  J  i  1  1  J 


(5) 


is  minimized.  The  least  squares  normal  equations  are  obtained  by  differentiating 
(5)  with  respect  to  bj  .  The  least  squares  normal  equations  are 


l  wj.k)  B(tj)  BT(t.)  b(k)* 
j*l  J  J  J 


l  ^k)  B(tj)  R(tj) 

j 


where  B^(tj)  is  the  vector  of  cubic  B-splines 
BT(t.)  *  [Bi(tj)  B 2 ( t j )  -  -  -  Bn ( tj ) 3 


(6) 


(7) 


Due  to  the  nature  of  the  B-splines  the  positive  definite  matrix  on  the  left  of 

(6)  is  banded  with  three  bands  above  and  below  the  main  diagonal.  To  conserve 

storage  the  four  distinct  diagonals  of  this  matrix  are  stored  as  columns  of  a 

''f  kl 

vertical  matrix.  The  dimension  of  the  vector  bv  '  is  HOTS  +  2  where  NOTS  is 
the  number  of  interior  knots.  The  banded  least  squared  normal  equations  are 
solved  by  a  banded  Cholesky  decomposition  algorithm.  The  sums  of  both  sides  of 
(6  )  are  performed  sequentially  so  that  all  of  the  ranges  and  weights  are  not 
needed  in  core  simultaneously.  The  IRWLS  can  be  continued  for  a  fixed  number 
of  iterations  or  until  the  fit  has  converged. 


INITIAL  WEIGHTS 

In  many  situations  the  IRWLS  procedure  works  successfully  when  all  of  the 
initial  weights  are  set  to  one,  i.e.,  the  iteration  is  started  with  an  ordinary 


7 


unweighted  least  sqwes  solution.  We  have  found  that  the  use  of  the  unweight¬ 
ed  least  squares  start  will  usually  result  in  convergence  of  the  IRWLS  cubic 
spline  to  a  good  fit  with  outliers  correctly  identified,  but  that  many  fewer 
iterations  are  required  if  a  robust  choice  of  initial  weights  is  used.  When 
outliers  are  present  in  either  the  first  or  last  intervals,  the  choice  of 
Initial  weights  in  these  intervals  is  most  important. 

The  initial  weights  for  the  robust  cubic  spline  fit  are  chosen  on  a  localized 
basis.  Let  R(t^),  1*1,  fiPTS(K)  be  the  range  measurements  in  the  K— knot 
interval.  To  determine  the  weights  wj°\  1=1,  NPTS(K)  in  the  K— interval  a 
linear  curve  is  robustly  fitted  to  the  measurements  in  the  interval.  Several 
methods  for  robustly  fitting  the  linear  curve  have  been  tried,  including  the 
nested  median  method  of  Siegel  [  1  ],  the  method  of  Theil  [  2],  a  modified  Theil 
method,  and  an  M-estimate  using  a  Hampel  ^-function.  Host  methods  performed 
about  equally  well  on  the  data  sequences  tested.  The  results  of  some  of  these 
tests  are  given  in  Appendix  A.  Because  of  its  simplicity,  the  modified  method 
of  Theil  was  selected  for  routine  application.  This  method  is  described  in  the 
following  paragraph. 

•  4L 

Let  R  be  the  median  of  the  observations  in  the  knot  interval. 

R  ■  median  (R(t. )} 

1-1,  NPTS1 (K)  (8 

Let  t  be  the  time  corresponding  to  R.  The  median  R  can  be  represented  as  the 
average  of  two  observations, 

5  •  <*<V)  *  R'W/Z  <! 

where  mi  »  m2  if  NPTS(K)  is  odd.  Define  the  set  of  slopes  {S.} 


c  ,  R(tJ  -  R 

a1  — - - — 

ti  -  t 


1*1,  NPTS(K) 
i  *  ma,  m2 


(10) 


Let  I  be  the  median  of  the  slopes, 

5  *  median  {S^} 

1  -  1 ,  NPTS(K) 

1  **  mj ,  m2 


Let  r.  be  the  residual , 


r . 

i 


R(^) 


5(ti  -  t),  1 


1,  NPTS(K) 
Let  r  be  the  median  of  these  residuals, 


01) 


(12) 


r  *  median  {r^} 
i  »  i,  NPTS(K) 

Now  compute  the  residuals. 


ri 


?1  -  r,  1  -  1,  NPTS(K) 
.(•) 


03) 

04) 


The  Initial  weights,  W}0',  are  computed  from  these  residuals  using  a  Hampel 
^-function. 


*C±) 

w(  o)  ,  Jv. 1  »  1. 

'  ft) 


NPTS(K) 


05) 


Where  s^  is  the  robust  dispersion  parameter, 

sk  »  median  (|ri|}/.6745  (16) 

Since  the  main  concern  in  setting  the  initial  weights  is  to  protect  the  cubic 
spline  fit  from  the  gross  outliers,  the  break  points  of  the  Hampel  ^-function 
are  set  at  a  »  2,  b  »  3,  c  »  4. 
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SOME  EXAMPLES 

Several  hundred  data  passes  have  been  run  with  the  fitting  procedure  described. 
Since  there  Is  an  average  of  maybe  five  receivers  on  each  data  pass,  the 
robust  preprocessing  method  described  has  been  used  on  more  than  one  thousand 
measurement  sequences.  The  method  has  performed  successfully  on  all  of 
these  sequences.  Most  of  these  sequences  are  rather  uneventful,  having  only 
a  few  Isolated  outliers.  There  have  been  some  sequences  which  have  some 
rather  dense  bursts  of  outliers.  These  sequences  best  Illustrate  the  ability 
of  the  method  described  to  detect  outliers.  Fig  1  presents  a  range  measurement 
sequence  and  Fig  2  the  robust  cubic  spline  fit  to  this  sequence.  Note  that  the 

outliers  In  Fig  1,  which  have  been  darkened,  occur  in  many  sizes.  The 
outliers  at  the  top  of  the  graph  were  added  by  hand  since  they  all  occurred 
far  off  scale  at  the  top.  The  sequence  of  Fig  1  has  about  10%  outliers. 

All  outliers  have  been  successfully  detected  and  removed  by  the  robust  spline 
fit.  The  measurements  in  Fig  1  have  two  dense  burst  of  outliers,  one  In  the 
Interval  (62356.6,  62359.5)  and  another  In  the  interval  (62367.2,  62375.4). 

The  measurement  sequence  in  Fig  3  has  outlier  bursts  In  the  intervals  (62358.4, 
62362.6),  (62374.4,. 62376.4),  and  (62379.7,  62382.4).  The  sequence  In  Fig  3 
has  about  15%  outliers.  The  sequence  In  Fig  5  has  bursts  of  outliers 
during  the  Intervals  (63117.8,  63122.6)  and  (63128.2,  63131.5).  Any  points 
away  from  the  main  curve  should  be  considered  outliers  In  Figs  1,  3,  and  5. 

Note  also  In  Figs  1,  3,  and  5  that  there  are  time  breaks  in  the  measurement 
sequences,  another  Important  consideration  In  preprocessing.  The  cubic 
spline  fit  to  the  sequence  of  Fig  1  Is  given  In  Fig  2.  The  cubic  spline 
fit  to  the  sequence  of  Fig  3  Is  given  In  Fig  4  and  the  cubic  spline  fit 
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to  the  measurement  sequence  In  Fig  5  Is  given  in  Fig  6.  The  knot  intervals 
in  Figs  2,  4,  6  are  designed  to  contain  twenty  data  points.  Note  that  some 
of  the  time  breaks  have  been  filled  with  fitted  data  points.  The  filling 
of  the  time  breaks  is  controlled  by  the  length  of  the  time  breaks  in 
relation  to  the  sample  rate  and  the  proportion  of  outliers  found  in  a 
knot  interval.  The  robust  cubic  spline  preprocessor  has  deleted  all 
outliers  from  the  measurement  sequence,  generated  measurements  during  the 
time  breaks  as  desired,  and  synchronized  different  measurement  sequences 
if  desired.  In  addition  the  measurement  variances  are  available  for 
further  processing.  The  IRWLS  cubic  spline  fit  converged  in  3  -  4 
iterations  for  the  examples  displayed.  This  fairly  rapid  convergence  is 
dependent  on  a  robust  method  for  choosing  good  Initial  weights.  Surprisingly, 
the  IRWLS  cubic  spline  iteration  for  these  examples  also  converges  using  an 
unweighted  least  squares  start,  but  at  the  expense  of  more  Iterations. 

For  the  measurement  sequences  displayed  here  the  IRWLS  cubic  spline  fit 
converged  in  7  -  8  iterations  using  an  unweighted  least  squares  start. 

Thus,  at  least  in  these  examples,  a  good  choice  of  the  initial  weights 
results  only  in  a  significant  improvement  In  computing  efficiency  and  not  in 
an  Improvement  of  fit.  Besides  a  good  selection  of  initial  weights,  another 
Important  choice  Is  the  number  of  data  points  per  knot  Interval,  NPTO. 

NPTO  must  be  large  enough  so  that  Is  likely  that  only  a  fraction,  say  less 

than  one  fourth  of  the  data  points  In  any  interval  will  be  outliers.  On 

the  other  hand,  If  NPTO  Is  too  large,  the  robust  linear  curve  fit  may  not 

be  a  good  enough  representation  of  the  variation  of  the  data  in  the  Interval. 
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This  appendix  describes  several  methods  of  choosing  the  Initial  weights 
for  the  robust  cubic  spline  preprocessing  and  compares  the  results  of  using 
these  methods  on  several  data  sets.  Each  of  these  methods  robustly  fits  a 
linear  curve  In  each  of  the  knot  Intervals  and  then  computes  the  initial 
weights  from  the  curve  fit  residuals  using  a  Hampel  ^-function.  Let 
R(tj),  t^  ■  1 ,  NPTS(k)  be  the  range  measurements  In  the  k— knot  Interval. 

Then  Method 


Define  the  slopes  s 


1J 


‘ij 


RCtj)  -  R(t^) 


V  *1 


j>1 


Let  i  be  the  median  of  these  slopes. 


s  »  median  {s,.} 
1.J  J 
J>1 


I 

—  Define  the  residuals  r^, 

ri  -  R(tf)  -  s  t, 

i 

Let  r  be  the  median  of  the  residuals,  r^ 

?  »  median  (r . } 

1»1 ,NPTS(k)  1 

Then  the  residuals  r^  ■  r^  -  r  are  used  to  compute  the  initial  weights  with  a 
Hampel  ^-function. 
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Nested  Medians 


Nested  or  repeated  medians  is  a  robust  regression  method  recently  described 
by  Siegel  [1].  Siegel  shows  that  this  method  has  the  highest  breakdown 
method  of  any  known  method.  This  method  Is  particularly  easy  to  apply  for 
a  linear  fit.  It  is  similar  to  the  Thell  method  and  modified  Thell 
method  already  described. 

Define  the  slopes  s^, 

R(t .)  -  R(t.) 

»u  •  Jtj  ■- 1,  <A-'> 


Define  1^  by 


51 


median  {s.,} 
j-1,NPTS(k)  1J 
j-1 


(A-2) 


and  further  let  s  be  defined  by 


5  •  median  {i^> 

1»1 ,NPTS(k)  1 


Similarly,  let  a^  be  the  Intercepts 


Mj 


RCt^tj  -  R(tjtj 


j-1 


Define  a^  as 


a.  ■  median  {a,.,} 

1  j-1,  NPTS(k) 

j-1 


( A— 3 ) 


(A-4) 


(A-S) 
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and  further  define  a  by 


I  ■  median  {a,}  ( A-6) 

1»1 ,NPTS(k)  1 

Let  r^  be  the  residuals 

ri  -  R(t.)  -  a  -  s  tr  1*1,  NPTS(k)  (A-7) 

The  weights  are  computed  from  these  residuals  using  a  Hampel 
function. 


The  following  data  sets  were  taken  from  the  knot  Intervals  of  the  data 
sequences  used  previously  to  Illustrate  the  application  of  the  robust 
range  measurement  preprocessing.  The  first  data  set,  shown  In  Fig  A1  is 
taken  from  the  measurement  sequence  given  In  Fig  1.  The  measurements  In  this 
set  are  from  the  time  Interval  62356.6  -  62359.5.  The  second  data  set, 
shown  In  Fig  A2  Is  taken  from  the  measurement  sequence  in  Fig  5.  The 
measurements  In  this  set  are  from  the  time  Interval  63128.2  -  63131.6. 

In  each  of  the  data  sets  the  weights  are  calculated  from  the  residuals  r^  by 


u(°) 

"l 


) 


) 


(4-8) 


where  *(•)  Is  a  Hampel  ^-function  with  breakpoints  2.,  3.,  4.  In  both  of 
those  data  sets  there  are  eight  outliers  In  the  sample  of  twenty.  Each  of 
the  robust  linear  methods  seem  to  have  no  difficulty  In  Identifying  the 
outliers  In  these  data  sets. 


FIGURE  A1 
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